###################################################################
#####check AFT mean estimates for individual simulation############
###################################################################
######################Figure S16###################################
k=1
numobs=100
setwd("c:/DASEV/S2/")

fname <- paste("simresult_N_",numobs,"_",k,".Rdata",sep="")
fname1 <- paste("otherresult_N_",numobs,"_",k,".Rdata",sep="")
fnamep <- paste("paradata_N_",numobs,"_",k,".Rdata",sep="")
fname.data <- paste("simdata_N_",numobs,"_",k,".Rdata",sep="")

load(file=fname.data)
datatmp <- data2
count <- rowSums(datatmp==0)
rawpall  <- count/(2*numobs)

load(file=fname) #sim.result
est <- sim.result$estimates
p <-  sim.result$pvalue
p.m <-  sim.result$pvalue_mean
post.mucontrol.f <- est[,1]
post.mucase.f <- est[,1]+est[,2]
post.pcontrol.f <- 1 - exp(est[,3])/(1+exp(est[,3]))
post.pcase.f <- 1 - exp(est[,3]+ est[,4])/(1+exp(est[,3])+est[,4])
post.var <- est[,5]^2

load(file=fname1) #other.result
t.p <- other.result[,1]
t.var <- other.result[,2]
t.mucontrol <- other.result[,3]
t.mucase <- other.result[,4]
t.pcontrol <- other.result[,5]
t.pcase <- other.result[,6]
t.p.m <- other.result[,7]
t.p.p <- other.result[,13]

aft.p <- other.result[,19]
aft.mucontrol <- other.result[,20]
aft.mucase <- other.result[,21]

twoT.mucontrol <- other.result[,28]
twoT.mucase <- other.result[,29]

load(file=fnamep) #parameters


sim.lod <- paradata[,1]
sim.mucontrol <- paradata[,2]
sim.pcontrol <- paradata[,3]
sim.sd <- paradata[,4]
sim.mucase <- paradata[,5]
sim.pcase <- paradata[,6]
sim.geneid <- paradata[,7]
sim.geneid1 <- paradata[,8]
numint <- paradata[,9]

sim.var <- sim.sd^2

q<-p.adjust(p,method='BH')
qmean<-p.adjust(p.m,method='BH')

t.q<-p.adjust(t.p,method='BH')
t.qmean<-p.adjust(t.p.m,method="BH")
par(mfrow=c(2,4))
par(oma=c(2,2,0,0))
par(mar=c(5,3,2,1))
ranges <-c(-4,4)
lfc.DA.TLK<- t.mucontrol[sim.geneid==1]-t.mucase[sim.geneid==1] 
hist(lfc.DA.TLK, breaks = 100,main="TLK DA", xlim = ranges,xlab = "")
lfc.DA.DASEV <- post.mucontrol.f[sim.geneid==1]-post.mucase.f[sim.geneid==1] 

hist(lfc.DA.DASEV, breaks = 100,main = 'DASEV DA', xlim = ranges,xlab = "")
lfc.DA.AFT <- aft.mucontrol[sim.geneid==1]-aft.mucase[sim.geneid==1] 
hist(lfc.DA.AFT, breaks = 100,main="AFT DA", xlim = ranges,xlab = "")
lfc.DA.twoT <- twoT.mucontrol[sim.geneid==1]-twoT.mucase[sim.geneid==1] 
hist(lfc.DA.twoT, breaks = 50,main="2T DA", xlim = ranges,xlab = "")

lfc.DA.TLK<- t.mucontrol[sim.geneid==0]-t.mucase[sim.geneid==0] 
hist(lfc.DA.TLK, breaks = 200,main="TLK non-DA", xlim = ranges,xlab = "")
lfc.DA.DASEV <- post.mucontrol.f[sim.geneid==0]-post.mucase.f[sim.geneid==0] 
hist(lfc.DA.DASEV, breaks = 100,main = 'DASEV non-DA', xlim = ranges,xlab = "")
lfc.DA.AFT <- aft.mucontrol[sim.geneid==0]-aft.mucase[sim.geneid==0] 
hist(lfc.DA.AFT, breaks = 100,main="AFT non-DA", xlim = ranges,xlab = "")
lfc.DA.twoT <- twoT.mucontrol[sim.geneid==0]-twoT.mucase[sim.geneid==0] 
hist(lfc.DA.twoT, breaks = 50,main="2T non-DA", xlim = ranges,xlab = "")

mtext("Log Fold Change", side=1, line=0, cex=2,outer=TRUE)
###################################################################